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NEW  HIGHER-ORDER  BOUNDARY-LAYER  EQUATIONS 
FOR  LAMINAR  AND  TURBULENT  FLOW  PAST  AXISYMMETRIC  BODIES 

C.  KLEINSTREUER*,  A.  EGHLIMA**  AND  J.E.  FLAHERTY'*’ 

Rensselaer  Polytechnic  Institute 
Troy,  NY  12181 

Abstract.  New  sets  of  boundary-layer  equations  accounting  for  flow  field 
non-uniformities  such  as  curvature  effects,  normal  stress  and  pressure  variations 
as  well  as  separation,  are  derived.  The  boundary-layer  flow  domain  is  subdivided 
into  (1)  a  parabolic  region  where  the  fluid  flow  is  approximately  parallel  to  the 
submerged  body,  i.e.  v«u  and  (2)  an  elliptic  region  which  includes  the  line  of 
separation  where  significant  interactions  between  the  boundary-layer  and  the  outer 
potential  flow  occur,  i.e.  v  »  u.  Closure  for  the  turbulent  flow  equations  has  to  be 
obtained  with  submodels  for  the  Reynolds  stresses  which  reflect  the  effects  of 
boundary-layer  thickening  as  well  as  separation.  The  accuracy  of  the  parabolic 
equations  was  compared  with  Van  Dyke's  higher-order  boundary- layer  equations  for 
laminar  flow  past  a  body  with  longitudinal  curvature.  The  results  demonstrate 
that  the  new  modeling  equations  make  a  measurable  difference  as  expected  from 
observations  made  by  Bradshaw  and  others. 

1.  Introduction.  An  accurate  but  also  computationally  efficient 
description  of  internal  or  external,  laminar  and  turbulent  flow  fields  is 
important  for  the  optimal  design  of  a  variety  of  mechanical  systems.  For 
example,  simulation  of  the  proposed  boundary-layer  equations  can  aid  in 
improving  the  design  of  bent  diffusors  and  submerged  bodies  in  terms  of 
reduced  pressure  loss  and  total  drag  reduction  respectively.  Furthermore, 
propeller  performance  of  marine  crafts  could  be  enhanced  with  a  more 
accurate  prediction  of  the  velocity  field  in  the  near-wake  region. 

At  the  high  Reynolds  numbers  of  interest,  the  analysis  involves  the 

usual  difficulties  encountered  with  turbulent  boundary  layers  in  a  pressure 
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gradient  now  accentuated  by  the  boundary  layer  thickness  (relative  to  the 
transverse  and/or  longitudinal  body  radius)  and  eventually  the  effects  of 
flow  separation.  Van  Dyke  (1962,  1969)  has  shown,  that  transverse  curvature, 
makes  a  contribution  to  differential  and  integral  flow  parameters  that  is 
additive  to  that  of  longitudinal  curvature  and  of  the  same  relative  order. 
Meroney  and  Bradshaw  (1975)  have  measured  turbulent  boundary  layer  growth 
in  a  prolonged  bend  and  show  that  a  small  (one  percent)  change  in  the  curva¬ 
ture  of  a  convex  surface  produces  a  relatively  large  (ten  percent)  change 
in  the  integral  properties  of  the  flow  field,  such  as  total  drag.  Further 
measurements  of  curvature  effects  on  laminar  and  turbulent  boundary-layer 
flow  parameters  were  published  by  Huang  et  al.  (1980),  Gillis  and  Johnson 
(1980),  Smith  et  al.  (1979),  Patel  and  Lee  (1978),  Ramaprian  and  Shivaprasd 
(1978),  and  Patel  (1974).  Different  investigators  (e.g.  Rastogi  and  Whitelaw 
(1971),  Bradshaw  (1973,  1975),  Patel  and  Lee  (1978),  Granville  (1978), 

Cebeci  et  al .  (1978),  Huang  et  al.  (1979),  Cebeci  (1979),  and  Patel  and 
Choi  (1980))  have  used  or  described  different  methods  to  solve  thick  bound¬ 
ary  layer  problems  for  special  situations  by  taking  into  account  some  of 
the  non-uniformities,  such  as  longitudinal  curvature  or  transverse  curvature 
effects,  but  not  all  of  them  as  discussed  earlier.  On  the  other  hand  we 
will  show  in  subsequent  sections  that  the  equation  usually  employed  for 
calculating  the  longitudinal  curvature  effect  (Van  Dyke  (1969))  is  not 
complete.  Our  derivations  (Sections  3  &  5)  indicate  that  additional  terms 
of  the  same  order  of  magnitude  as  Van  Dyke's  second  order  terms  are  necessary 
for  representing  the  influence  of  longitudinal  curvature  on  the  growth  of 
boundary  layers.  Authors  (e.g.  Cebeci  (1971),  Rastogi  and  Whitelaw  (1971), 
Cebeci  et  al .  (1978,  (1979))  who  employed  a  turbulent  version  of  Van  Dyke's 
equation  could  not  predict  the  curvature  effects  on  integral  properties  of 


the  flow  field  as  measured  by  Meroney  and  Bradshaw  (1975).  Solutions  of  the 
boundary  layer  equations  employing  prescribed  variations  of  the  pressure 
gradient  (Patel  (1974))  failed  near  the  point  of  flow  separation  and  con¬ 
sequently  a  singular  behavior  in  the  boundary  layer  solution  was  postulated 
(Goldstein  (1948)).  However,  Williams  (1977)  used  the  results  of  various 
researchers  (e.g.  Keller  and  Takami  (1966),  Son  and  Hanratty  (1969),  Dennis 
(1970),  and  Masliyah  (1970)),  who  solved  the  Navier-Stokes  equation  for  flows 
that  include  the  point  of  zero  stress,  to  show  that  such  a  singular  behavior 
is  not  a  physical  property  of  the  flow  but  it  is  a  characteristic  of  the 
solutions  of  the  boundary  layer  equations. 

We  are  presenting  higher-order  boundary-layer  equations  that  simulate 
incompressible  laminar  and  turbulent  flow  fields  more  accurately  than  with 
existing  models. 

2.  System  conceptualization  and  approach.  Consider  the  steady  two- 
dimensional  or  axisymmetric  flow  of  a  Newtonian  fluid  past  a  stationary  body. 
Figure  1  schematically  depicts  the  various  flow  developments  including  the 
three  major  regions  of  interest  at  the  tail  of  the  submerged  body  where  the 
boundary  layer  is  thick  due  to  curvature  effects.  In  region  I  the  mean  flow 
streamlines  remain  nearly  parallel  to  the  solid  surface  regardless  of  the 
relative  thickness  of  the  boundary  layer.  In  region  II  streamlines  are  not 
parallel  to  the  body  surface  and  the  velocity  component  u  parallel  to  the 
body  surface  is  of  the  same  order  of  magnitude  as  the  velocity  component  v 
normal  to  the  body.  The  flow  in  this  region  is  characterized  by  strong 
interactions  between  the  boundary  layer  and  the  outer  inviscid  flow  in  region 
III.  Basically  two  new  sets  of  governing  equations  for  thick  turbulent 


Flow  Fields  of  Interest 


boundary  layer  flows  (regionl  and  II)  are  derived.  Laminar  flow  fields  with 
any  type  of  curvature  effect  can  be  regarded  as  a  special  case  of  the  first 
complete  set  of  modeling  equations. 

The  approach  for  simulating  a  general  external  or  internal  flow  system  is 

1 

best  illustrated  with  reference  to  Fig.  1.  Starting  from  the  stagnation  point, 
a  suitable  computer  code  (e.g.  Cebeci  et  al .  (1978)) is  employed  until  the  boundary 
layer  becomes  thick,  i.e.  «  1  is  not  valid  any  more.  The  intermediate  results 
are  then  relayed  to  the  program  of  our  new  parabolic  equations  for  region  I  which 
are  accurate  as  long  as  v  «  u  which  may  include  the  point  of  separation.  We 
note  that  boundary  layer  thickening  might  occur  for  certain  body  geometries 
at  quite  an  early  stage  as  demonstrated  in  Section  6  for. a  laminar  flow  case 
study.  When  reverse  flow  occurs  and  v  =  u,  the  new  elliptic  model ing  equa¬ 
tions  of  region  II  are  activated.  Interaction  with  the  potential  flow 
(Cebeci  et  al .  (1978))  is  simulated  with  the  local  displacement  thickness 
using  an  iterative  procedure.  In  order  to  preserve  computational  efficiency, 
a  continuous  function  is  assumed  for  the  pressure  drop  from  the  calculated 
pressure  at  the  separation  point  to  the  wake  pressure  behind  the  afterbody 
which  is  apparently  of  the  order  of  magnitude  of  the  outer  flow  pressure. 

The  near-wake  velocity  field  is  then  actually  computed  with  Schl ichting’s  solu¬ 
tion  for  (circular)  wakes  behind  a  three-dimensional  body  (e.g.  White  (1974), 

Riley  and  Metcalfe  (1980)).  More  detailed  investigations  for  wake  flow  are  dis¬ 
cussed  by  Wu  (1972),  Pope  and  Whitelaw  (1976)  and  recently  by  Smith  (1979). 

3.  Derivation  of  the  governing  equations  for  region  I.  The  governing 
equations  for  the  generalized  class  of  thick  boundary  layer  developments  are 
written  in  terms  of  fluid  mass  and  momentum  fluxes  (cf.  Kleinstreuer  (1982)): 


(3.1) 


(3.2) 


||  +  7-pv  =  0 


at 


P  V  +  V*pv  V  =  -  V-TT  +  pf 


where  pv  is  the  mass  flux  (i.e.,  fluid  mass  per  unit  area  and  per  unit  time), 
pvv  is  the  momentum  flux  and  ir  =  p6  +  t  is  the  total  stress  tensor  consisting 
of  the  thermodynamic  pressure  as  well  as  shear  plus  normal  stresses. 

We  restrict  our  analysis  to  the  steady  two-dimensional  or  axi  symmetric 
flow  of  a  Newtonian  fluid  around  a  submerged  body  and  introduce  the  tangential 
and  normal  (s,  n)  orthogonal  coordinate  system  shown  in  Fig.  2.  Under  these 
conditions,  the  equations  of  continuity  and  momentum  may  be  written  as: 


(3.3) 
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where  the  geometric  parameters  h  =  1  +  jj-  and  K  =  i  reflect  longitudinal 
curvature  effects  whereas  r  =  r0  + n  cos8  accommodates  the  transverse  curvature 
(cf.  Figure  2). 


h  =  1  +  Kn  . 


The  system  of  equations  (3. 3-3. 5)  can  be  simplified  for  distinct  parts  of  the 
flow  domain  employing  the  relative  order  of  magnitude  analysis. 

In  region  I  of  Fig.  2,  the  normal  component  of  the  velocity  vector  v  is 
significantly  smaller  than  the  longitudinal  component  u.  With  this  assumption, 
we  introduce  dimensionless  variables  that  are  of  the  order  of  unity  in  the 
boundary  layer  (cf.  Van  Dyke  (1969))  as 

u  =  u/U  ,  v  =  Rem  v/U 

(3.6)  s  =  s/L  ,  n  =  Rem  n/L 

Re  =  UL/\>  ,  p  =  p/p  U2 

Non-dimensional ization  of  the  governing  equations  is  achieved  by  inserting 
the  dimensionless  variables  (3.6)  into  first  the  continuity  equation  (3.3). 
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This  equation  after  restoring  original  variables  and  rearrangement  will 
reduce  to 


(3.9) 
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Substituting  (3.6)  into  equation  (3.4)  yields: 

U2  u  3u  ,  U2  -  3D  +  Ui  kT  Re~m 
L  l  +  KnRe'm3§  L  3R  L  1  +  i<n  Re_m  09  = 

_ui _ J _ 3|+  v{jy__L__i!u  + 

L  l  +  KnRe~m  3s  L2  1 +  ICn  Re“m  3s2  . 


(3.10) 


+  M  +  JL  iL^.  £  _ 1  _ 5r 

L2  af1  L2  3s  L(r0+n  Re"mcos6)n  + j<n  Re_Tn)2  3s 

- __! - %  ]  +  XL  Re2m  9u  r  — _ - _J _ 

( 1  +  kTn  Re-rn) 3  gn  L2  Sn  L(r0+n  Re~ni  cose)3n 


+  _ — ,, — i - ah.  ]•> 

(l  +  K  n  Re_m)  3n  '  . 


Using  the  relation 
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one  can  find  for  the  previous  equation 
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After  some  rearrangements,  the  equation  of  motion  in  s-direction  can  be 
written  as 


(3.12) 
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Restoring  the  original  variables  in  equation  (3.12)  will  yield 
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Substituting  (3.6)  into  equation  (3.5)  yields 
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Using  again  the  relation  yyy  =  1-c  +  0(e2)  one  obtains 


(3.15) 
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After  rearrangements  and  simplifications,  we  have 
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Restoring  the  original  variables  yields 

(3.17)  Ku*  =  1|£+  o(Re_1,  Re’2”,  •••)  . 

p  on 


Hence,  for  region  I,  (3.9,  3.13  and  3.17)  are  the  basic  equations  where  the 
flow  indices  i  and  j  indicate  the  cases  i=j  =  0  for  flow  past  a  flat  plate  and 
i  =  l,  j  =  o  for  axisymmetric  bodies  with  no  longitudinal  curvature  and  i=0, 
j  =  l  for  flow  without  transverse  curvature.  Finally,  i  =  j  =  t  when  longitudinal 
and  transverse  curvature  effect  exist.  The  exponent  m  is  known  for  turbulent 
thin  boundary  layer  flow  with  zero  pressure  gradient  along  a  flat  plate  (m  = 

1  9m 

see  Schl ichting,  1979,  p.  638).  Thus,  terms  of  the  order  of  Re"  ,  Re  ,  etc. 


are  neglected. 

Depending  on  the  geometry  of  the  submerged  body  it  is  possible  that  the 
boundary-layer  thickness  60  is  not  small  (i.e.  =  0(1)  )  but  the  mean  flow 

streamlines  still  remain  nearly  parallel  to  the  surface  so  that  the  normal 
component  of  the  velocity  v  is  still  much  smaller  than  the  longitudinal 
component  u.  In  this  case,  dimensionless  variables  that  are  of  the  order  of 
unity  in  the  boundary  layer  are  introduced  as: 
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Inserting  the  dimensionless  variables  into  the  system  (3. 3-3. 5),  and  following 
the  same  procedure  that  was  used  previously,  one  obtains 


4?  (u  r)  +  Re-m  (vh  r)  =  0 
(3.19)  a||  +  R  Re"m  v|S+uvK=-||  +  OfRe"1) 

Re'm  ti  ||  +  *  J!  =  -  ||  +  0(Re-2m.  Re'"1-1)  . 

Restoring  the  original  variables 

£(ur)  +  ^  (vh  r)  =  0 

<3'20)  “f+"vf,uv^-If+0(Re-1l. 

u||-Ku2=--1|^  +  0(Re“2m,  Re~m_1) 


For  the  turbulent  flow  analysis,  time-averaged  variables  are  denoted  by 
overbars  and  random  fluctuations  by  primes,  so  that 

(3.21)  u=u+u'  ,v=v+v'  and  p  =  p  +  p' 

where  Reynolds  stresses  generated  by  the  non-axi symmetric  velocity  fluctuations 
w'  are  neglected. 

In  addition,  the  following  rules  of  operation  for  time-averaged  variables 
are  required  for  reference 

f  =  f ,  TT^  =  f  +  g ,  TTq  *  f-g,  —  =  ||  and  7fds  =  /fds 
Hence,  u2  =  u2  +  u'2  and  u*v  =  u*v  +  u ‘v1 


Equations  (3.22  to  3.24)  contain  five  unknowns,  namely,  u,  v,  p,  u* v' ,  (\p s ) . 

In  order  to  gain  closure  it  is  therefore  necessary  to  furnish  some  additional 
relationships. 

Based  on  the  objective  of  computational  efficiency,  and  with  the  measured 
data  sets  found  in  the  open  literature,  we  postulate  an  algebraic  turbulence 
or  zero-equation  model  for  the  shear  stresses  reflecting  curvature  effects  and 
relate  these  expressions  directly  to  the  normal  stresses.  Among  the  algebraic 
submodels,  the  mixing  length  hypothesis  (MLH)  has  proven  very  useful  and 
realistic  for  a  wide  range  of  case  studies,  provided  that  good  choices  were 
made  for  the  mixing  length  distributions.  However,  the  MLH  cannot  represent, 
in  general,  recirculating  flows,  turbulence  convection  or  diffusion  and  non-zero 
effective  viscosities  or  Prandtl  numbers  at  the  points  of  zero  velocity  gradients. 
Two  commonly  employed  transport  equations  for  turbulent  shear  stresses  are  the 
mixing  length  formula  (Von  KSrm5n  (1931) ) 


(3.27a) 


-u 1 v 1  =  Z 


m 


3u 


3u 

3y 


and  the  eddy  viscosity  formula  (Boussinesq  Q877)) 


(3.27b) 


3u_ 
't  dy 


where  =  vt(£m)  again. 

Following  the  method  used  for  thin  turbulent  boundary  layer  calculations, 
a  composite  layer  consisting  of  two  regions  is  postulated  (e.g.  Baker  &  Launder 
0974),  Kwon  &  Pletcher  (1979)).  The  distributions  of  £m  and  vt  are  then 
described  by  two  separate  empirical  expressions.  The  criterion  to  define 
the  inner  and  outer  region  of  the  boundary  layer  is  the  continuity  of  the 
eddy  viscosity  o£  the  mixing  length  at  the  hypothetical  "interface".  The 


shear  stresses  are  a  function  of  mean-time  velocity  and  velocity  gradient  as 
well  as  mixing  lengths  which  in  trun  are  dependent  upon  the  radii  of 
curvatures,  viz: 

(3.28)  uV  =  uV  [vt(im,R,r0) ;  u,Vu]. 

The  turbulent  flow  analysis  as  well  as  predicted  vs.  measured  results  can  be 
found  in  Eghlima  (1982)  and  will  be  published  in  a  forthcoming  paper. 


4.  The  governing  equations  for  region  II.  In  region  II  which  possibly 
includes  the  point  of  separation,  we  can  no  longer  assume  that  v  «  u.  Hence,, 
the  flow  in  this  region  is  characterized  by  strong  interaction  between  the 
(detached)  boundary  layer  and  the  outer  potential  flow.  It  is  also  conceiv¬ 
able  that  the  attached  (parabolic)  boundary  layer  regime  may  be  influenced  by 
the  (elliptic)  reverse  flow  region  after  separation.  Actual  simulation  results 
of  thick  separating  boundary  layer  flows,  viscid-inviscid  interactions  and 
near-wake  effects  will  be  discussed  in  a  future  series  of  papers.  In  this 
section  we  summarize  the  governing  equations  for  region  II. 

After  the  appropriate  correction  in  equation  (3.6),  since  now  v  *  u, 
and  substitution  into  the  system  of  equations  (3.3  to  3.5)  we  obtain 

(4.1)  ^  u)  +  Mr’-r’O-lr*))]  +  0(Re~2m)  =  0 


(4.2) 


uM+v|£+uvK  =  _l|£+vliu  + 

oj  on  p  dS 


v  +  0(Re"2m) 
9n 


(«.3) 


We  note  that  transverse  curvature  effects  only  appear  in  the  continuity 
equation  (4.1);  whereas  the  approximated  s-  and  n-momentum  equations  (4.2 
and  4.3)  only  explicitly  contain  the  longitudinal  curvature  parameter  in 
one  term.  However,  in  contrast  to  the  first-order  boundary  layer  equations, 
the  Navier-Stokes  equation  in  n-direction  is  almost  fully  preserved  since 
v  =  u.  In  addition,  transverse  curvature  effects  have  to  be  coupled  with 
the  Reynolds  stresses. 

Again  splitting  the  instantaneous  variables  into  time-smoothed  and 
random  components  we  can  derive  the  steady  turbulent  flow  equations  for 
region  II  as 

,(4.4)  £  (rj  2)  *  £  (M  v)  =0 


(4.5) 


(4.6) 


where  M  =  rJ  -  rQ(l-hJ)  and  K  =  (^)^J"  . 

The  associated  boundary  conditions  and  the  meaning  of  i  and  j  are  given 


—  3u  .  —  3u  .  — 
u— — +v  —  +  uvK  « 
ds  on 


-±!*  +  v 

p  3s 


32u  1  ,  ,M  +  roN  3 


3n 


a  [(±  2 


)  -r-  («■  )  + 

05  • 


+  sin  0(uf  )  +  -^  (u'v'  M)  +  MK  (u'v')] 


n|5+7||-Kj2  .-i|S  +  vp.i{^(u.v-ro)  + 


3s  ‘  3n 


+  2)  (T72)  +  (cos  0)  Cv’2)  +  K  r0[(v*2)  -  (u'2))} 


in  Section  3. 


>.  Case  study  for  laminar  boundary-layer  flow  with  longitudinal 
curvature  effect.  To  demonstrate  the  technical  merit  of  the  new  governing 
equations,  a  spinoff  from  equation  (3.4),  the  special  case  of  laminar  boundary 
layer  flow  with  longitudinal  curvature  effects  is  analyzed.  Van  Dyke  (196Z), 
0.969)  developed  a  higher-order  boundary  layer  equation  which  takes  into 
account  the  effect  of  longitudinal  curvature  as 


(5.1) 


u|i+hjv|i+Kuv,.l|£  +  vhj2ii  ,vit|a. 


3n 


p  3s 


3n“ 


3n 


In  contrast,  our  s-momentum  equation  which  can  be  directly  obtained  from  (3.13)  for 
reads 


u  It  +  hJ  v  f£+  Kuv  +  u  f£  +  (1-hJ)  v  |^  ]  = 


3u 


3u 


3u 


(5.2) 


3S 


3n 
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3n 


'  -  ?  i 4  v  hj  fr 4  v  K !? 4  i  If 4  vo-if  >  0  ]• 


Hence,  the  proposed  modeling  equation  (5.2)  carries  four  additional  terms 
which  are  of  the  same  (second)  order  Of  magnitude  when  compared  with  equation 
(5.1).  For  h=l ,  which  includes  K=0,  both  equations  reduce  to  Prandtl’s 
boundary  layer  equation. 

Before  both  equations  are  numerically  compared  it  is  of  interest  to  trace 
the  derivation  procedures  leading  to  the  noted  discrepancy.  The  traditional 

departure  point  for  obtaining  a  higher-order  approximation  of  laminar  incom¬ 
pressible  flow  along  a  curved  surface  is  either  the  system  of  equations  (7) 
on  page  119  in  Goldstein  (1948)  or  equations  (22  and  23)  on  page  202  in 
Rosenhead  (1963).  The  two  sets  of  equations  are  mathematically  equivalent 
but  physical  significance  is  lost viien  Goldstein's  equation  is  multiplied  by 
H  (original  notation)  to  generate  Rosenhead's  equation.  Indeed,  Rosenhead's 


equation  can  be  obtained  in  multiplying  equation  (3.4)  by  h  =  1  +  and  then 
neglecting  terms  of  the  order  of  (■£■)*  and  higher.  During  this  procedure 
certu.n  terms  of  the  order  of  ( ^ )  could  not  be  retained.  This  is  best 
illustrated  by  considering  an  equation  of  the  functional  form 

(5.3a)  ^  +  Be  +  0(e2)  + - =0 

where  h  =  1  +  e  +  0(e2)  and  e  =  ^  .  • 

Equation  (5.3a)  is  representative  for  Goldstein's  equation  whereas  Rosenhead's 
equation  would  take  on  the  form 


(5.4a) 

A  +  Bhe  +  0(e2)  +  —  =  0 

or 

(5.4b) 

A  +  Be  +  fl(e2)  +  ...  =  0. 

This  last  procedure,  leading  to  Van  Dyke's  equation  (5.1),  indicates  a  short¬ 
coming  when  compared  with  equation  (5.3)  now  rewritten  with  h~^  =  (1+e)”^  = 
1-e  +  0(cz ) 

(5.3b)  (1-e)  A  +  Be  +  0(e2)  +  •••  =  0. 

In  equation  (5.3b)  which  is  representative  of  our  equation  (5.2),  the  term, 

-e  A ,  is  preserved  in  the  derivation  procedure. 

In  Section  6,  numerical  results  from  a  case  study  document  the  differ¬ 
ences  between  equations  (5.2)  and  (5.1). 


6.  Preparation  of  governing  equations  and  numerical  comparison.  In 
order  to  simplify  the  comparison  of  numerical  case  studies,  equations  (3.23) 
and  (3.24)  are  rewritten  in  a  general  form. 


(6.1) 


A,  u  ^  +  A2  v  -1^-  +  A-,  uv  =  A4  +  v  Ac  (A,  f^)+vA.  —■ 
1  3s  on  J  ^  3s  ->  Sn  o  ?n  7  3n 


(6.2) 


A3  u 


2  =  A  & 
h8  Sn 


A  — 
m9  3n 


Now,  Van  Dyke's  equation  (5.1)  is  obtained  by  setting 

A1  =  1  ,  A2  =  hj  ,  A3  =  K  ,  A4  =  -  1  ,  A5  =  h^  ,  Ag  =  1  and  A?  =  K  . 

Our  s-momentum  equation  (5.2)  requires 

A]  =  2-  ,  A2  =  1  ,  A3  =  K  ,  A4  =  -(2-h)/p  ,  A5  =  Ag  =  1  and  A?  =  K  . 


A  finite  difference  method  developed  by  H.B.  Keller  and  described  by 
Cebeci  and  Smith  (1974)  was  selected  after  improvements  were  implemented  as  the 
best  candidate  for  our  modeling  equations  (Eghlima  0982)).  In  order  to  prepare 
the  governing  equations,  submodels  and  boundary  conditions  for  the  numerical 
solution  procedure,  the  following  steps  have  to  be  accomplished:  (1)  trans¬ 
formation  of  the  modeling  equations  into  a  new  coordinate  system  (modified 
Mangler-Levy-Lee  transformation),  (2)  transformation  of  the  momentum  equations 
and  associated  boundary  conditions  into  a  system  of  first-order  partial  differ¬ 
ential  equation,  and  (3)  discretization  of  the  modeling  equations  into  finite 
difference  forms  using  centered-differenced  quotients.  Some  salient  aspects  of 
the  first  two  steps  are  given  below. 

In  order  to  remove  the  large  variations  of  the  boundary  layer  thickness 
and  to  obtain  the  governing  equations  (3.22  to  3.24)  in  simple  form,  an  appro¬ 
priate  coordinate  transformation  is  implemented. 

(6.3a)  d  (  =  Pe  p  t  ue  ds 

where,  in  general,  the  turbulence  viscosity  =  ye(l  +  em)  and  =  em(5,n); 

dn  =  [pe  ue/  /2£~  ]  (  £“  dn 

r 


(6.3b) 


In  using  the  continuity  equation  (3.22),  the  stream  function  can  be 
defined  as 


(6.4a,b) 


u  = 


rV 


3n 


and 


1  3£ 

“  Mpe  3s 


A  dimensionless  stream  function  f(£,n)  is  related  to  ,n)  as: 
(6-5)  iMs,n)  =  /ITl/'fU.n)  . 


Now  the  partial-derivative  operators 
For  example: 


can  be  calculated. 


(6.6a)  u  =  u^f'  ,  where  f '  = 

*  e  9rj 

(6.6b)  v  =  -  ve  ue  /TTL^/H  t  £  +  ||  ♦  (ff  )n  f’]  . 


Substituting  the  transformed  variables  and  their  derivatives  into  equations 
(3.23  and  3.24)  -  equation  (3.22)  is  automatically  satisfied  -  the  momentum 
equations  in  the  g,  n )-  coordinate  system  read: 


(6.7) 


[Q1  f"]'  +  Q2  f’f"  +  Q3(f')2  +  Q4  ff"  +  Q5  ff'  +  Qfi  f"  = 


+  Qc 


f" 


f' 


3f 

9£ 


(6.8) 


Q„  f"  +  Q12(f’)2  -  Q13(ff )? 


where  the  Q’s  are  functions  of  geometric,  fluid,  and  flow  parameters.  In 
substituting  (  )^  from  equation  (6.8)  into  (6.7),  a  single  equation  results. 

To  solve  this  third  order,  nonlinear  partial  differential  equation  with 
Keller's  Box  Method,  new  dependent  variables  U(£,n)  and  V(£,n)  are  introduced 
and  a  coupled  system  of  first-order  PDE's  is  generated: 


(6.9a) 

(6.9b) 


(6.9c) 


f'  =  U  =  u/ue 
U‘  =  V 


[Q-J  V]'  +  q2  uv  +  Q14(U)2  +  q4  f  v  +  q5  f  u  +  Q15  V  = 


(ff) 

^  T) 


+  Q8  u  ||  Mo  V 


3i 

as 


0  u  — 

W10  u 


The  system  (6.9)  constitutes  the  new  set  of  higher-order  boundary  layer  equations 
for  laminar  and  turbulent  flows  in  region  I. 

The  associated  boundary  conditions,  formerly  given  as  (3.25),  can  be  written  as: 
(6.10a)  f’(£,0),  i.e.  U(£,0)  =  0  and  •  f(S_,0)  =  0. 


In  addition, 
(6.10b) 

where 


n  -»■  ,  f' 


1 


1  +  D(s,n) 


D(C.n)  =  k 


cos  e 


For  the  case  of  flow  past  a  flat  plate,  D(5,n)  =  0  or  f'  =1.  The  edge 
boundary  condition  may  be  sensitive  to  the  specification  of  n^.  To  reduce 
this  sensitivity,  the  expression  (6.10b)  is  differentiated  so  that  with 


(6.10c) 

where 


n  -  0m  ,  f"  =  - 


0,(5, n) 


[1  +D(€,n)]' 


2Lr  COS  0  STr~  -1/2 


In  combining  (6.10b  and  c),  the  outer  boundary-layer  condition  reads 

(6. lOd)  n  -  nro  ,  f"  +  Dj (c,n)  f’2  =  0 

V  +  D,U2 


or 


=  0. 


The  coefficients  of  (6.9c)  are  defined  as 
*14 
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.  n  n  1l2  3n 

Q3  -  q7  5^  35 
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n  n  5,1  °n 
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-  25  Q5 

'Ve^ 

^12 

-  a3/7T 

h  • 


Q13  "  ‘  A8  Ug  ^  Lr  ^  and  A7  =  G(s’n)  =  ^  y  3n  +  ei/a  M  s1nS  +  M  3n  + 


For  the  simulation  of  (5.1)  and  (5.2),  G(s,n)  =  K  where  r0-*°°  and  e_  h  0,  and 
the  appropriate  parameters  A-|  to  Ag  have  to  be  inserted. 

The  boundary  conditions  for  (5.1)  and  (5.2)  are 

U(5,n  =  0)  «  0  =  f(c,  n  =  0)  and  U(5,n-nJ~  1. 

Following  the  procedure  for  implementing  Keller’s  Box  Method  as  outlined  in 
Cebeci  and  Smith  (1974),  equations  (6.9a-c)  are  programmed  and  then  solved 
with  an  efficient,  generalized  equation  solver  developed  by  Varah  (1976)  and 
implemented  by  Flaherty  and  Mathon  (1980).  The  pressure  gradient  is  obtained 
from  the  solution  of  the  inviscid  flow  regime  (region  III)  for  which  the 
code  developed  by  Cebeci  et  al.  (1978)  was  employed. 


To  demonstrate  the  effect  of  longitudinal  curvature  and  to  show  the  numerical 
differences  between  equations  (5.1)  and  (5.2),  we  selected  an  axisymmetric 
body  designed  and  documented  by  Huang  et  al .  (1979).  Figure  3  depicts  schem¬ 
atically  one-half  of  the  symmetric  body  which  in  our  case  is  semi-infinite 
since  r0  -*■  Of  much  higher  resolution  and  physical  interest  is  Fig.  4.  It 
shows  the  local  ^  =  K 6  vs.  the  chord  length  of  the  submerged  body.  The  curve 
peaks  early  due  to  a  sharp  decrease  of  the  radius  of  longitudinal  curvature  and 

then  falls  to  zero  (despite  the  growing  boundary  layer)  when  R  When  the 

* 

flat  middle  section  is  over,  increases  again  since  R  becomes  finite  and  6 

is  quite  thick  (=  0.4  ft)  by  then.  The  threshold  value  for  ( -| )  =  0.003 

K  o 

was  suggested  by  Cebeci  and  Smith  (1974)  based  on  Bradshaw's  measurements  which, 

JP  T  —  * 

however,  were  made  for  turbulent  flows.  For  values  ^  >  (^)  curvature 
effects  become  significant. 

Figure  5  demonstrates  that  our  "additional  second-order  terms"  have  indeed  a 
significant  effect,  actually  10%  at  one  point  in  this  case  study.  It  is 

6 

evident  that  the  maximal  difference  occurs  after  (x)  has  peaked.  There  are 
two  reasons  for  this  shift.  Most  influential  is  the  term  containing  the  pressure 
gradient  which  is  retained  in  our  derivation.  In  addition,  there  is  an  "after¬ 
effect"  of  the  peak  on  the  velocity  distribution  (enhanced  by  our  four  additional 
terms)  which  becomes  fully  visible  a  small  distance  downstream.  The  influential 
pressure  gradient  is  about  equal  to  zero  at  the  maximum  point  for  (CFA-CFV)/CFF 
which  explains  the  sudden  drop  of  the  curve  somewhat  in  correspondence  with  ^  (x). 
Tov/ards  the  tail  of  the  body  becomes  significant  again  and,  indeed, 
(CFA-CFV)/CFF  f  0. 


o 
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Fig.  3.  Semi-infinite  body  with  longitudinal  curvature  (adapted  from  Huang  et  al .  1980) 
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Fig.  4.  Local  =  K6  v$.  chord  length  of  submerged  body. 
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Fig.  5.  Effects  of  thick  laminar  boundary  layer  development  on  the  local  friction 

coefficient:  A  comparison  between  equations  (5.1)  and  (5.2)  at  Re  =  (?(10-s). 


In  contrast  to  Fig.  5  where  modeling  improvements  for  points  along  the  solid 

surface  are  shown,  i.e.  „  (x).  Figs.  6a-c  depict  the  differences  in  boundary 

dy  wall 

layer  velocity  profiles  as  produced  by  (5.1)  and  (5.2).  Again,  the  higher 
accuracy  of  our  equation  is  quite  evident  and  has  a  significant  effect  on  the 
local  fluid  flow  parameters,  especially  near  and  at  the  wall, as  well  as  on  the 
point  of  separation.  Figure  7  shows  the  development  of  the  pressure  gradient 
as  discussed  earlier. 

7.  Conclusions  and  future  work.  Based  on  a  mathematically  rigorous 
analysis,  new  sets  of  higher-order  boundary  layer  equations  for  axisymmetric 
laminar  and  turbulent  flow  fields  were  derived.  For  the  case  of  thick  laminar 
boundary  layer  flow  with  longitudinal  curvature  effects,  Van'Dyke's  momentum 
equation  was  numerically  compared  with  our  new  equation  which  carries  four  addi¬ 
tional  (second-order)  terms.  The  difference  between  the  two  equations  in  terms 
of  local  as  well  as  integral  properties  of  the  flow  system  is  significant.  It 
is  anticipated  that  the  existing  discrepancy  betv/een  predicted  and  measured 
data  for  thick  turbulent  flows  can  now  be  closed  with  the  use  of  the  new 
higher-order  boundary  layer  equations.  Hence,  several  case  studies  will  be 
conducted  which  will  include  the  simulation  of  laminar  and  turbulent  boundary 
layer  developments  along  longitudinally  and  transversely  curved  surfaces  and 
comparison  with  measured  data  sets  published  by  So  and  Mellor  (1973),  Simpson 
et  al.  (1974)  and  Huang  et  al .  (1979,  1980). 

The  turbulent  flow  data  collected  by  Simpson  et  al.  (1974)  include  separation. 

It  is  obvious  that  separating  flows  are  very  difficult  to  investigate.  Strong 
interactions  between  the  potential  flow  and  thick  turbulent  boundary  layer 
have  to  be  accounted  for.  In  addition,  separation  and  near-wake  effects  might 
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influence  the  upstream  region  of  the  still  attached  boundary  layer.  However, 
we  will  test  the  performance  of  our  parabolic  equations  of  region  I  also  for 
regions  containing  mild  recirculation  so  that  the  complex  simulation  of  (4.4 
to  4.6)  can  be  avoided  whenever  possible. 

In  any  case,  only  simple  submodels  for  the  Reynolds  stresses  will  be  employed 
in  order  to  demonstrate  also  for  turbulent  flow  the  physical  significance  of 
our  new  equations.  Previously,  empirical  turbulence  models  for  specific  case 
studies  were  used  in  conjunction  with  incomplete  momentum  equations  to  match 
laboratory  observations. 

To  reduce  the  cumbersome  and  error-prone  work  of  system  discretization,  the 
mesh  generator  developed  by  Kleinstreuer  (1980)  and  advancements  proposed  by 
Dwyer  et  al .  (1980)  will  be  implemented  for  future  case  studies. 
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NOMENCLATURE 


to  Ag  parameters  in  equations (6. 1  &  6.2) 

f  dimensionless  stream  function  defined  by  Eq.  (6.5) 

h.|Eh=l  +  Kn  metric  coefficient  (dimensionless) 
h2  =  1  metric  coefficient  (dimensionless) 

h3  =  r  metric  coefficient  (dimensionless) 

K  =  1  longitudinal  curvature  (ft)-1 

Lr  reference  length  (ft) 

mixing  length  (ft) 

M=  (r)1  -  (ro)1  [l-(h)J]  (ft) 
n  coordinate,  normal  to  body  (ft) 

P  pressure 

p  non-dimensional  pressure 

Q1  to  Q15  parametersin  equation  (6.9c) 

R  radius  of  longitudinal  curvature  (ft) 

r =  r0 + n  cos0  transverse  radius  of  curvature  (ft) 
r0  transverse  radius  of  curvature  of  the  body  (ft) 

s  coordinate  along  the  body  surface  (ft) 

s  non-dimensional  coordinate  along  the  body  surface 


t  time 

u  £  component  of  the  velocity  (ft/sec) 

u  dimensionless  velocity  in  .s  direction 

u  time-averaged  velocity  in  s_ direction  (ft/sec) 

u'  velocity  fluctuation  in  s_  direction  (ft/sec) 

u'2  Reynolds  normal  stress  in  direction 


dimensionless  variable 


uh|I 

3n 

Uinf 


v 

v 

V 


v'2 

V  =  2U 
v  “  an 


w 


w 


u‘v 


free  stream  velocity 
£  component  of  the  velocity  (ft/sec) 
dimensionless  velocity  in  n_  direction 
time-averaged  velocity  in  n  direction  (ft/sec) 
velocity  fluctuation  in  n_  direction  (ft/sec) 
Reynolds  normal  stress  in  ji  direction 
non-dimensional  variable 
velocity  normal  to  (s-n)  plane  (ft/sec) 
Reynolds  normal  stress  normal  to  (s^n)  plane 
Reynolds  shear  stress 


Greek  Symbols 

3  longitudinal  curvature  effect  parameter 

6  angle  between  axis  of  symmetry  and  tangent  to  the 

surface  (radiant) 

v  kinematic  viscosity 

Vj.  eddy  viscosity 

V  dynamic  viscosity  (lbm/ft-sec) 

P  mass  density  (lbm/ft3) 

K  transformed  s-coordinate  (lbm/ft-sec) 


transformed  n-coordinate  (non- 
boundary-layer  thickness  (ft) 
non-dimensional  eddy  viscosity 


Superscripts 

i  transverse  curvature  index 

=1  does  exist 
=0  does  not  exist 

j  longitudinal  curvature  index 

=1  does  exist 
=0  does  not  exist 


Subscripts 


e 

eff 

f 

t 


outer  edge  of  boundary  layer 
effective  value 
flat  plate 
turbulent 


1 


dimensional ) 


w 


wall 


